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1  Activities  and  results  reported  in  the  pre¬ 
vious  report  -  October  2006 

Two  series  of  investigation  were  reported  in  the  progress  report  of  October 

2006. 

The  first  one  was  the  study  of  spherical  shock  waves  (blast  waves)  prop¬ 
agating  in  dusty  gases  and  flow  behind  them.  The  purpose  of  the  study 
was  to  see  if  we  would  observe  the  phenomena  in  which  particles  drive  the 
gas.  At  the  beginning  when  a  spherically  expanding  shock  wave  is  strong, 
solid  particles  that  are  originally  stationary  are  accelerated  by  the  gas  behind 
the  shock  wave.  Later  when  the  shock  wave  expands  further,  it  decelerates 
and  becomes  weaker.  At  this  stage  the  particles  have  already  gained  certain 
amount  of  momentum  from  the  gas  behind  the  shock  wave.  Consequent  ly, 
under  certain  conditions,  it  can  be  expected  that  particles  now  accelerate  the 
gas,  contrarily  to  the  early  stage  of  the  explosion. 

In  conclusion,  however,  we  did  not  observe  the  aforementioned  phenom¬ 
ena  at  least  within  the  range  of  parameters  used  in  the  current  numerical 
analyses.  The  conclusion  at  the  time  of  the  previous  report  was  that  the 
range  of  parameters  such  as  initial  shock  strength,  particle  diameters  etc., 
was  not  wide  enough.  It  was  felt  that  more  simulations  were  necesary  to 
make  any  decisive  conclusion. 

In  the  second  series  of  investigations,  the  non-equilibrium  region  behind 
a  plane  shock  wave  is  studied  by  numerically  integrating  the  conservation 
equations  of  mass,  momentum  and  energy.  Unsteady  drag  forces  obtained 
both  experimentally  and  numerically  were  modeled  and  implemented  in  the 
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numerical  code.  It  is  reported  that  the  drag  force  is  about  an  order  of  mag¬ 
nitude  larger  compared  to  the  steady  one  when  a  shock  wave  passes  over  a 
sphere  due  to  complex  wave  interactions.  The  purpose  of  the  investigation 
was  to  evaluate  the  effect  of  unsteady  drag  force  on  the  structure  of  non¬ 
equilibrium  region  behind  a  plane  shock  wave.  The  unsteady  drag  forces 
were  obtained  for  stationary  solid  spheres  of  different  diameters  (Sun  et  al. 
2004).  In  a  dusty  gas,  a  particle  starts  to  move  at  the  moment  when  a  shock 
wave  arrives  at  the  particle.  Simulations  had  been  carried  out  for  both  cases 
of  with  and  without  the  particle  movement  and  the  results  were  compared. 
This  series  of  computations  had  not  been  completed  in  the  previous  report. 
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2  Progress  report  on  the  activities  after  the 
previous  report 

Following  activities  have  been  carried  out  since  the  previous  progress  report. 

1.  Investigations  on  the  effect  of  unsteady  drag  force  on  the  non-equilibrium 
region  behind  plane  shock  wave. 

2.  Development  of  three  dimensional  numerical  code  for  simulating  flows 
in  a  dusty  gas. 

3.  Investigations  of  spherically  expanding  shock  waves  in  dusty  gases. 

The  item  1  is  the  continuation  of  the  research  activities  described  in  the 
previous  section.  Computations  of  drag  force  for  different  incident  shock 
Mach  numbers  were  done  to  investigate  the  effect  of  incident  shock  Mach 
number  on  the  unsteady  drag  force.  The  effect  of  different  mass  ratio  between 
gas  and  particle  phases  was  also  investigated. 

The  results  were  presented  at  the  26th  international  shock  wave  sympo¬ 
sium  (” Effect  of  unsteady  drag  force  on  structure  of  non-equilibrium  region 
behind  shock  wave  in  gas-particle  mixture”  by  T.  Saito,  M.  Saba,  M.  Sun 
and  Iv.  Takayama,  Proc.  26th  ISSW  ).  Also  a  paper  entitled  ’’The  effect  of 
an  unsteady  drag  force  on  the  structure  of  a  non-equilibrium  region  behind 
a  shock  wave  in  a  gas-particle  mixture”  is  submitted  to  the  international 
journal  of  Shock  Waves  and  is  accepted  for  publication. 

The  manuscripts  of  the  proceedings  of  the  ISSW26  and  the  manuscript 
submitted  to  the  Shock  Wave  Journal  are  included  in  this  progress  report  as. 
Both  will  be  published  soon. 
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The  second  item  has  practical  purpose  of  applications  to  numerical  sim¬ 
ulations  of  several  environmental  problems  such  as  volcano  eruptions,  dust 
explosions  etc.  with  appropriate  terrains  into  considerations.  Investigations 
of  generation  and  propagation  of  blast  waves  by  explosive-type  volcano  erup¬ 
tions  have  been  carried  out  in  our  group  for  many  years.  We  expect  that 
the  development  of  the  3D  numerical  code  will  provide  us  more  realistic  in¬ 
formation  to  study  the  mechanism  of  volcano  eruptions  and  also  to  establish 
better  safety  measures. 

The  third  item  is  a  continuation  of  the  work  that  is  described  in  the 
previous  section.  Due  to  limited  range  of  computational  parameter,  we  did 
not  observe  the  phenomenon  that  particles  drive  the  gas.  However,  we  felt 
that  there  is  much  to  be  studied  on  behavior  of  blast  waves  in  dusty  gases. 
Simulations  with  wider  range  of  computational  parameters  are  now  being- 
considered  and  simulations  are  being  carried  out. 
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3  Development  of  3D  numerical  code 


Three-dimensional  numerical  code  is  developed  for  simulating  shock  wave 
propagation  in  dusty  gases.  The  code  just  started  to  run  and  we  are  testing 
the  code  by  applying  it  to  a  couple  of  simple  test  problems. 

The  code  is  briefly  described  in  this  section  by  showing  the  basic  equations 
that  are  solved  numerically.  Numerical  results  of  test  problems  are  also 
presented. 

3.1  Basic  equations  and  numerical  method 

The  governing  equations  of  the  dusty  gas  flows  in  three  spatial  dimensions 
(3D)  can  be  written  as: 
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and  for  solid  particles 
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The  vector  of  interaction  terms  is 
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E  =  p  jc^T  +  |  (?r  +  V2  +  w 2)  |  ,  (6) 

Vt  =  a  |cm0  +  i  ( ud 2  +  vd2  +  wd)  | ,  (7) 

and  here  u  and  ud  are  the  x-component  of  gas  and  dust  velocity  and  similarly 
v  and  vd  for  the  v-component.  The  x  and  y-component  of  the  drag  force  are 
given  as 

D x  =  \nd2p  {u  -  ud)  | u  -  ud |  CD,  (8) 

O 

Dv  =  ^cl2p  ( v  -  vd)  |t  -  vd\  CD.  (9) 

Dz  =  \nd2p(w  -  wd)  | w  -  wd |  CD.  (10) 

O 

Here  as  in  the  ID  case,  the  following  models  of  the  drag  coefficient  and  the 
Nusselt  number  are  assumed. 


CD  =  0.48  +  28Re-0'85, 

Nu  =  2.0  +  0.6Pi'3Re2, 

The  particel  Reynolds  number  is  defined  in  2D  as  follows, 

Re  _  pd\j {u  -  ud)2  +  (v  -  vd )2  +  (w  -  wd )2 

/' 
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As  to  the  solution  strategy,  the  operator  splitting  is  used  for  handling  the 
interaction  terms.  The  homogeneous  part  of  the  governing  equations  (14) 
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and  (15)  below  are  solved  first  and  then  the  solution  is  further  modified  by 
taking  the  effect  of  interaction  between  the  gas  and  solid  particle  phases  by 
solving  (16)  and  (17). 
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In  solving  (14)  and  (15),  we  had  options  to  use  Strang  type  operator  splitting 
to  treat  the  multi  space  dimensions  (Strang,  G.  (1968):  On  the  construction 
and  comparison  of  difference  schemes,  SIAM  J.  Numer.  Anal.  5(3):506-517). 
In  this  study,  however,  the  numerical  fluxes  in  both  space  dimensions  are 
considered  at  once. 

For  the  evaluation  of  numerical  fluxes  of  the  particle  phase,  flow  param¬ 
eters  are  interpolated  between  neighboring  computational  cells  as  follows, 
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here  q°  is  the  initial  value  at  the  cell  interface, 
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The  interpolated  values  of  flow  parameters  at  the  cell  interfaces  are  ob- 
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tained  from  the  foilwing  equations 
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The  numerical  fluxes  at  cell  interfaces  are  calculated  from  the  values  with 
t  =  |A t.  Full  details  of  the  computational  methods  for  2D  are  available 
in  the  reference  (Saito,  T.,  Marumoto,  M.,  Takayama,  K.(2003):  Numerical 
investigations  of  shock  waves  in  gas-particle  mixtures.  Shock  Waves,  No. 13, 
299-322). 

3.2  Preliminary  test  calculations 

Two  problems  were  calculated  for  testing  the  numerical  code.  Although,  in 
principle,  the  code  is  simply  an  extension  of  our  2D  code,  it  is  significantly 
more  complex.  The  code  is  developed  on  general  curvilinear  coordinates  and 
runs  in  parallel  on  multiple  processing  elements  by  using  MPI  libraries. 

3.2.1  Explosion  of  a  pressurezed  cube  under  a  duty  gas  layer 

The  initial  state  of  the  problem  is  presented  in  Fig.  1.  A  cubic  container 
filled  with  high  pressure  gas  without  dust  particles  is  placed  on  the  ground. 
The  region  between  the  two  gray  horizontal  planes  is  filled  with  a  dusty  gas. 
The  container  is  broken  at  time  zero  generating  a  blast  wave.  The  blast  wave 
interacts  with  the  dustv-gas  layer  and  the  state  of  the  mass  concentration  at 
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some  time  later  is  displayed  in  Fig.  2.  It  is  seen  that  the  borders  between 
pure  gas  and  the  dusty  gas  are  deformed  by  the  blast  wave  and  the  flow 
behind  it. 

Since  this  is  only  a  test  for  checking  the  code  performance,  not  much 
to  discuss  from  the  physical  point  of  view.  However,  the  three  dimensional 
nature  of  the  flow  is  confirmed.  Also  we  can  see  that  the  parallel  execution 
of  the  program  is  performed  nicely,  since  there  is  no  obvious  disturbances 
observed  near  the  border  of  the  regions  calculated  with  different  processors. 


Figure  1:  Initial  state  of  the  test  problem,  r  =  0.0. 


3.2.2  Explosion  of  a  high  pressure  dusty  gas  in  a  cubic  container 

In  contrast  to  the  previous  problem,  a  high  pressure  dusty  gas  is  exploded 
into  pure  gas.  Figure  3  presents  the  initial  condition  of  the  problem. 

After  the  instantaneous  burst  of  the  cubic  container  that  is  filled  with 
a  high  pressure  dusty  gas  at  some  distance  above  the  ground,  blast  wave  is 
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Figure  2:  Distribution  of  particle  mass  concentration,  r  =  20.0. 

generated  and  propagates  into  the  surrounding  pure  gas.  The  distribution 
of  particle  mass  concentration  at  a  normalized  time  (r  =  20.0)  is  shown  in 
Fig.  4.  Complex  distribution  pattern  of  solid  particles  from  the  container  is 
observed. 

In  conclusion,  it  seems  that  the  3D  numerical  code  is  developed  and 
running  properly. 
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Figure  3:  Initial  state  of  the  explosion  of  a  cubic  containar  filled  with  dusty 
gas,  t  =  0.0. 


Figure  4:  Distribution  of  particle  mass  concentration  at  values  of  0.0001, 
0.001  and  0.2,  r  =  20.0. 


4  Future  work 


We  plan  to  carry  out  following  research  subjects. 

1.  Numerical  investigations  of  spherical  shock  waves  in  dusty  gases. 

2.  Applications  of  2D  and  3D  numerical  codes  to  practical  problems. 

As  a  continuing  research  project,  we  investigate  spherically  expanding 
shock  waves  in  dusty  gases.  We  anlvze  the  numerical  results  that  have  been 
accumulated  since  the  beginning  of  this  research  project. 

One  of  the  candidates  for  application  of  the  2D  code  is  the  simulation 
of  the  starting  and  stopping  of  supersonic  dusty  gas  flow  in  a  nozzle.  The 
subject  is  investigated  for  pure  gases  but  not  much  is  reported  for  cases  with 
dusty  gases. 

As  an  application  of  the  3D  numerical  code,  we  will  simulate  imaginary 
eruption  of  a  volcano.  With  the  3D  dusty  gas  code,  it  is  expected  that  we 
will  be  able  to  obtain  some  useful  information  regarding  pyrocrastic  flows. 
We  have  detailed  information  about  Mount  Usu,  one  of  active  volcanoes  in 
Hokkaido  Japan,  together  with  high  resolution  digital  elevation  map  of  the 
area. 
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